home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / SCSL / dtrsen.z / dtrsen
Encoding:
Text File  |  2002-10-03  |  11.7 KB  |  331 lines
  1.  
  2.  
  3.  
  4. DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DTRSEN - reorder the real Schur factorization of a real matrix A =
  10.      Q*T*Q**T, so that a selected cluster of eigenvalues appears in the
  11.      leading diagonal blocks of the upper quasi-triangular matrix T,
  12.  
  13. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  14.      SUBROUTINE DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S,
  15.                         SEP, WORK, LWORK, IWORK, LIWORK, INFO )
  16.  
  17.          CHARACTER      COMPQ, JOB
  18.  
  19.          INTEGER        INFO, LDQ, LDT, LIWORK, LWORK, M, N
  20.  
  21.          DOUBLE         PRECISION S, SEP
  22.  
  23.          LOGICAL        SELECT( * )
  24.  
  25.          INTEGER        IWORK( * )
  26.  
  27.          DOUBLE         PRECISION Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( *
  28.                         ), WR( * )
  29.  
  30. IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
  31.      These routines are part of the SCSL Scientific Library and can be loaded
  32.      using either the -lscs or the -lscs_mp option.  The -lscs_mp option
  33.      directs the linker to use the multi-processor version of the library.
  34.  
  35.      When linking to SCSL with -lscs or -lscs_mp, the default integer size is
  36.      4 bytes (32 bits). Another version of SCSL is available in which integers
  37.      are 8 bytes (64 bits).  This version allows the user access to larger
  38.      memory sizes and helps when porting legacy Cray codes.  It can be loaded
  39.      by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
  40.      only one of the two versions; 4-byte integer and 8-byte integer library
  41.      calls cannot be mixed.
  42.  
  43. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  44.      DTRSEN reorders the real Schur factorization of a real matrix A =
  45.      Q*T*Q**T, so that a selected cluster of eigenvalues appears in the
  46.      leading diagonal blocks of the upper quasi-triangular matrix T, and the
  47.      leading columns of Q form an orthonormal basis of the corresponding right
  48.      invariant subspace.
  49.  
  50.      Optionally the routine computes the reciprocal condition numbers of the
  51.      cluster of eigenvalues and/or the invariant subspace.
  52.  
  53.      T must be in Schur canonical form (as returned by DHSEQR), that is, block
  54.      upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2
  55.      diagonal block has its diagonal elemnts equal and its off-diagonal
  56.      elements of opposite sign.
  57.  
  58.  
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  71.  
  72.  
  73.  
  74. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  75.      JOB     (input) CHARACTER*1
  76.              Specifies whether condition numbers are required for the cluster
  77.              of eigenvalues (S) or the invariant subspace (SEP):
  78.              = 'N': none;
  79.              = 'E': for eigenvalues only (S);
  80.              = 'V': for invariant subspace only (SEP);
  81.              = 'B': for both eigenvalues and invariant subspace (S and SEP).
  82.  
  83.      COMPQ   (input) CHARACTER*1
  84.              = 'V': update the matrix Q of Schur vectors;
  85.              = 'N': do not update Q.
  86.  
  87.      SELECT  (input) LOGICAL array, dimension (N)
  88.              SELECT specifies the eigenvalues in the selected cluster. To
  89.              select a real eigenvalue w(j), SELECT(j) must be set to w(j) and
  90.              w(j+1), corresponding to a 2-by-2 diagonal block, either
  91.              SELECT(j) or SELECT(j+1) or both must be set to either both
  92.              included in the cluster or both excluded.
  93.  
  94.      N       (input) INTEGER
  95.              The order of the matrix T. N >= 0.
  96.  
  97.      T       (input/output) DOUBLE PRECISION array, dimension (LDT,N)
  98.              On entry, the upper quasi-triangular matrix T, in Schur canonical
  99.              form.  On exit, T is overwritten by the reordered matrix T, again
  100.              in Schur canonical form, with the selected eigenvalues in the
  101.              leading diagonal blocks.
  102.  
  103.      LDT     (input) INTEGER
  104.              The leading dimension of the array T. LDT >= max(1,N).
  105.  
  106.      Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
  107.              On entry, if COMPQ = 'V', the matrix Q of Schur vectors.  On
  108.              exit, if COMPQ = 'V', Q has been postmultiplied by the orthogonal
  109.              transformation matrix which reorders T; the leading M columns of
  110.              Q form an orthonormal basis for the specified invariant subspace.
  111.              If COMPQ = 'N', Q is not referenced.
  112.  
  113.      LDQ     (input) INTEGER
  114.              The leading dimension of the array Q.  LDQ >= 1; and if COMPQ =
  115.              'V', LDQ >= N.
  116.  
  117.      WR      (output) DOUBLE PRECISION array, dimension (N)
  118.              WI      (output) DOUBLE PRECISION array, dimension (N) The real
  119.              and imaginary parts, respectively, of the reordered eigenvalues
  120.              of T. The eigenvalues are stored in the same order as on the
  121.              diagonal of T, with WR(i) = T(i,i) and, if T(i:i+1,i:i+1) is a
  122.              2-by-2 diagonal block, WI(i) > 0 and WI(i+1) = -WI(i). Note that
  123.              if a complex eigenvalue is sufficiently ill-conditioned, then its
  124.              value may differ significantly from its value before reordering.
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  137.  
  138.  
  139.  
  140.      M       (output) INTEGER
  141.              The dimension of the specified invariant subspace.  0 < = M <= N.
  142.  
  143.      S       (output) DOUBLE PRECISION
  144.              If JOB = 'E' or 'B', S is a lower bound on the reciprocal
  145.              condition number for the selected cluster of eigenvalues.  S
  146.              cannot underestimate the true reciprocal condition number by more
  147.              than a factor of sqrt(N). If M = 0 or N, S = 1.  If JOB = 'N' or
  148.              'V', S is not referenced.
  149.  
  150.      SEP     (output) DOUBLE PRECISION
  151.              If JOB = 'V' or 'B', SEP is the estimated reciprocal condition
  152.              number of the specified invariant subspace. If M = 0 or N, SEP =
  153.              norm(T).  If JOB = 'N' or 'E', SEP is not referenced.
  154.  
  155.      WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  156.              On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  157.  
  158.      LWORK   (input) INTEGER
  159.              The dimension of the array WORK.  If JOB = 'N', LWORK >=
  160.              max(1,N); if JOB = 'E', LWORK >= M*(N-M); if JOB = 'V' or 'B',
  161.              LWORK >= 2*M*(N-M).
  162.  
  163.              If LWORK = -1, then a workspace query is assumed; the routine
  164.              only calculates the optimal size of the WORK array, returns this
  165.              value as the first entry of the WORK array, and no error message
  166.              related to LWORK is issued by XERBLA.
  167.  
  168.      IWORK   (workspace) INTEGER array, dimension (LIWORK)
  169.              IF JOB = 'N' or 'E', IWORK is not referenced.
  170.  
  171.      LIWORK  (input) INTEGER
  172.              The dimension of the array IWORK.  If JOB = 'N' or 'E', LIWORK >=
  173.              1; if JOB = 'V' or 'B', LIWORK >= M*(N-M).
  174.  
  175.              If LIWORK = -1, then a workspace query is assumed; the routine
  176.              only calculates the optimal size of the IWORK array, returns this
  177.              value as the first entry of the IWORK array, and no error message
  178.              related to LIWORK is issued by XERBLA.
  179.  
  180.      INFO    (output) INTEGER
  181.              = 0: successful exit
  182.              < 0: if INFO = -i, the i-th argument had an illegal value
  183.              = 1: reordering of T failed because some eigenvalues are too
  184.              close to separate (the problem is very ill-conditioned); T may
  185.              have been partially reordered, and WR and WI contain the
  186.              eigenvalues in the same order as in T; S and SEP (if requested)
  187.              are set to zero.
  188.  
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  203.  
  204.  
  205.  
  206. FURTHER DETAILS
  207.      DTRSEN first collects the selected eigenvalues by computing an orthogonal
  208.      transformation Z to move them to the top left corner of T.  In other
  209.      words, the selected eigenvalues are the eigenvalues of T11 in:
  210.  
  211.                    Z'*T*Z = ( T11 T12 ) n1
  212.                             (  0  T22 ) n2
  213.                                n1  n2
  214.  
  215.      where N = n1+n2 and Z' means the transpose of Z. The first n1 columns of
  216.      Z span the specified invariant subspace of T.
  217.  
  218.      If T has been obtained from the real Schur factorization of a matrix A =
  219.      Q*T*Q', then the reordered real Schur factorization of A is given by A =
  220.      (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the
  221.      corresponding invariant subspace of A.
  222.  
  223.      The reciprocal condition number of the average of the eigenvalues of T11
  224.      may be returned in S. S lies between 0 (very badly conditioned) and 1
  225.      (very well conditioned). It is computed as follows. First we compute R so
  226.      that
  227.  
  228.                             P = ( I  R ) n1
  229.                                 ( 0  0 ) n2
  230.                                   n1 n2
  231.  
  232.      is the projector on the invariant subspace associated with T11.  R is the
  233.      solution of the Sylvester equation:
  234.  
  235.                            T11*R - R*T22 = T12.
  236.  
  237.      Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote the
  238.      two-norm of M. Then S is computed as the lower bound
  239.  
  240.                          (1 + F-norm(R)**2)**(-1/2)
  241.  
  242.      on the reciprocal of 2-norm(P), the true reciprocal condition number.  S
  243.      cannot underestimate 1 / 2-norm(P) by more than a factor of sqrt(N).
  244.  
  245.      An approximate error bound for the computed average of the eigenvalues of
  246.      T11 is
  247.  
  248.                             EPS * norm(T) / S
  249.  
  250.      where EPS is the machine precision.
  251.  
  252.      The reciprocal condition number of the right invariant subspace spanned
  253.      by the first n1 columns of Z (or of Q*Z) is returned in SEP.  SEP is
  254.      defined as the separation of T11 and T22:
  255.  
  256.                         sep( T11, T22 ) = sigma-min( C )
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268. DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  269.  
  270.  
  271.  
  272.      where sigma-min(C) is the smallest singular value of the
  273.      n1*n2-by-n1*n2 matrix
  274.  
  275.         C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
  276.  
  277.      I(m) is an m by m identity matrix, and kprod denotes the Kronecker
  278.      product. We estimate sigma-min(C) by the reciprocal of an estimate of the
  279.      1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) cannot
  280.      differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
  281.  
  282.      When SEP is small, small changes in T can cause large changes in the
  283.      invariant subspace. An approximate bound on the maximum angular error in
  284.      the computed right invariant subspace is
  285.  
  286.                          EPS * norm(T) / SEP
  287.  
  288.  
  289. SSSSEEEEEEEE AAAALLLLSSSSOOOO
  290.      INTRO_LAPACK(3S), INTRO_SCSL(3S)
  291.  
  292.      This man page is available only online.
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  
  302.  
  303.  
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.                                                                         PPPPaaaaggggeeee 5555
  328.  
  329.  
  330.  
  331.